library(ashr)
library(cjoint)
library(tidyverse)
library(cregg)
library(xtable)

###SET WORKING DIRECTORY and LOAD FILES ----


set.seed(12345678)
df <- read.csv("b_data/3_long_data_for_analysis.csv")
white_long <- read.csv("b_data/4_white_data_for_analysis.csv")
black_long <- read.csv("b_data/5_black_data_for_analysis.csv")


####Reshaping data as factors ----

df$race <- as.factor(df$race)

df <- df %>% 
  mutate(Faces = factor(Faces, levels = c("White man", "White woman", "Light man", "Light woman", 
                                          "Medium man",   "Medium woman", "Brown man", 
                                          "Brown woman", "Black man", "Black woman"), 
                        ordered = is.ordered(Faces)))


df <- df %>%
  mutate(Background = factor(Background, levels = c("Legal immigrant",
                                                    "Undoc. immigrant",
                                                    "Naturalized citizen",
                                                    "US-born legal parents",
                                                    "US-born undoc. parents",
                                                    "US-born legal grandparents",
                                                    "US-born undoc. grandparents")))

df <- df %>%
  mutate(Identifies = factor(Identifies, levels = c("American",
                                                    "Latino",
                                                    "Hispanic",
                                                    "Mexican-American",
                                                    "Mexican",
                                                    "White",
                                                    "Black")))

df <- df %>%
  mutate(Language = factor(Language, levels = c("Only English",
                                                "Bilingual",
                                                "Limited English",
                                                "No English"),
                           ordered = is.ordered(Language)))

df <- df %>%
  mutate(Education = factor(Education, levels = c("High school dropout",
                                                  "High school graduate",
                                                  "Bachelor's degree",
                                                  "Advanced/graduate degree"),
                            ordered = is.ordered(Education)))

df <- df %>%
  mutate(Profession = factor(Profession, levels = c("Waiter", "Janitor",
                                                    "Sales manager", "Teacher",
                                                    "IT professional", "Doctor",
                                                    "Small business owner")))

df <- df %>%
  mutate(Religion = factor(Religion, levels = c("Presbyterian", "Evangelical",
                                                "Catholic", "Not religious")))

df <- df %>%
  mutate(Political = factor(Political, levels = c("Democrat", "Republican",
                                                  "Independent", "Not political")))

df <- df %>%
  mutate(Spouse = factor(Spouse, levels = c("Not married",
                                            "S. is American",
                                            "S. is Black",
                                            "S. is White",
                                            "S. is same as profile")))

black_long <- black_long %>% 
  mutate(Faces = factor(Faces, levels = c("White man", "White woman", "Light man", "Light woman", 
                                          "Medium man",   "Medium woman", "Brown man", 
                                          "Brown woman", "Black man", "Black woman"), 
                        ordered = is.ordered(Faces)))


black_long <- black_long %>%
  mutate(Background = factor(Background, levels = c("Legal immigrant",
                                                    "Undoc. immigrant",
                                                    "Naturalized citizen",
                                                    "US-born legal parents",
                                                    "US-born undoc. parents",
                                                    "US-born legal grandparents",
                                                    "US-born undoc. grandparents")))

black_long <- black_long %>%
  mutate(Identifies = factor(Identifies, levels = c("American",
                                                    "Latino",
                                                    "Hispanic",
                                                    "Mexican-American",
                                                    "Mexican",
                                                    "White",
                                                    "Black")))

black_long <- black_long %>%
  mutate(Language = factor(Language, levels = c("Only English",
                                                "Bilingual",
                                                "Limited English",
                                                "No English"),
                           ordered = is.ordered(Language)))

black_long <- black_long %>%
  mutate(Education = factor(Education, levels = c("High school dropout",
                                                  "High school graduate",
                                                  "Bachelor's degree",
                                                  "Advanced/graduate degree"),
                            ordered = is.ordered(Education)))

black_long <- black_long %>%
  mutate(Profession = factor(Profession, levels = c("Waiter", "Janitor",
                                                    "Sales manager", "Teacher",
                                                    "IT professional", "Doctor",
                                                    "Small business owner")))

black_long <- black_long %>%
  mutate(Religion = factor(Religion, levels = c("Presbyterian", "Evangelical",
                                                "Catholic", "Not religious")))

black_long <- black_long %>%
  mutate(Political = factor(Political, levels = c("Democrat", "Republican",
                                                  "Independent", "Not political")))

black_long <- black_long %>%
  mutate(Spouse = factor(Spouse, levels = c("Not married",
                                            "S. is American",
                                            "S. is Black",
                                            "S. is White",
                                            "S. is same as profile")))

black_long$amerimport_high <- as.factor(black_long$amerimport_high)

white_long <- white_long %>% 
  mutate(Faces = factor(Faces, levels = c("White man", "White woman", "Light man", "Light woman", 
                                          "Medium man",   "Medium woman", "Brown man", 
                                          "Brown woman", "Black man", "Black woman"), 
                        ordered = is.ordered(Faces)))


white_long <- white_long %>%
  mutate(Background = factor(Background, levels = c("Legal immigrant",
                                                    "Undoc. immigrant",
                                                    "Naturalized citizen",
                                                    "US-born legal parents",
                                                    "US-born undoc. parents",
                                                    "US-born legal grandparents",
                                                    "US-born undoc. grandparents")))

white_long <- white_long %>%
  mutate(Identifies = factor(Identifies, levels = c("American",
                                                    "Latino",
                                                    "Hispanic",
                                                    "Mexican-American",
                                                    "Mexican",
                                                    "White",
                                                    "Black")))

white_long <- white_long %>%
  mutate(Language = factor(Language, levels = c("Only English",
                                                "Bilingual",
                                                "Limited English",
                                                "No English"),
                           ordered = is.ordered(Language)))

white_long <- white_long %>%
  mutate(Education = factor(Education, levels = c("High school dropout",
                                                  "High school graduate",
                                                  "Bachelor's degree",
                                                  "Advanced/graduate degree"),
                            ordered = is.ordered(Education)))

white_long <- white_long %>%
  mutate(Profession = factor(Profession, levels = c("Waiter", "Janitor",
                                                    "Sales manager", "Teacher",
                                                    "IT professional", "Doctor",
                                                    "Small business owner")))

white_long <- white_long %>%
  mutate(Religion = factor(Religion, levels = c("Presbyterian", "Evangelical",
                                                "Catholic", "Not religious")))

white_long <- white_long %>%
  mutate(Political = factor(Political, levels = c("Democrat", "Republican",
                                                  "Independent", "Not political")))

white_long <- white_long %>%
  mutate(Spouse = factor(Spouse, levels = c("Not married",
                                            "S. is American",
                                            "S. is Black",
                                            "S. is White",
                                            "S. is same as profile")))

white_long$amerimport_high <- as.factor(white_long$amerimport_high)

black_long$raceimport_high <- as.factor(black_long$raceimport_high)




###CONJOINT DESIGN ----

attribute_list <- list()

attribute_list[["Faces"]] <-c("White man", "White woman", "Light man", "Light woman",
                              "Medium man", "Medium woman", "Brown man", "Brown woman",
                              "Black man", "Black woman")
attribute_list[["Background"]] <- c("Legal immigrant", "Undoc. immigrant",
                                    "Naturalized citizen", 
                                    "US-born legal parents",
                                    "US-born undoc. parents",
                                    "US-born legal grandparents",
                                    "US-born undoc. grandparents")
attribute_list[["Identifies"]] <- c("American", "Latino",
                                    "Hispanic", "Mexican-American",
                                    "Mexican", "White", "Black")
attribute_list[["Language"]] <- c("Only English",
                                  "Bilingual",
                                  "Limited English",
                                  "No English")
attribute_list[["Education"]] <- c("High school dropout",
                                   "High school graduate",
                                   "Bachelor's degree",
                                   "Advanced/graduate degree")
attribute_list[["Profession"]] <- c("Waiter", "Janitor",
                                    "Sales manager", "Teacher",
                                    "IT professional", "Doctor",
                                    "Small business owner")
attribute_list[["Religion"]] <- c("Presbyterian", "Evangelical",
                                  "Catholic", "Not religious")
attribute_list[["Political"]] <- c("Democrat", "Republican",
                                   "Independent", "Not political")
attribute_list[["Spouse"]] <- c("Not married", "S. is American",
                                "S. is Black", "S. is White",
                                "S. is same as profile")



constraint_list <- list()

constraint_list[[1]] <- list()
constraint_list[[1]][["Education"]] <- c("High school dropout",
                                         "High school graduate")
constraint_list[[1]][["Profession"]] <- c("Doctor",
                                          "IT professional", "Teacher")

constraint_list[[2]] <- list()
constraint_list[[2]][["Education"]] <- c("Bachelor's degree")
constraint_list[[2]][["Profession"]] <- c("Doctor")


constraint_list[[3]] <- list()
constraint_list[[3]][["Background"]] <- c("US-born legal parents",
                                          "US-born undoc. parents",
                                          "US-born legal grandparents",
                                          "US-born undoc. grandparents")
constraint_list[[3]][["Language"]] <- c("No English")

constraint_list[[4]] <- list()
constraint_list[[4]][["Background"]] <- c("Legal immigrant",
                                          "Naturalized citizen",
                                          "Undoc. immigrant")
constraint_list[[4]][["Language"]] <- c("Only English")

constraint_list[[5]] <- list()
constraint_list[[5]][["Profession"]] <- c("IT professional",
                                          "Doctor",
                                          "Sales manager",
                                          "Teacher")
constraint_list[[5]][["Language"]] <- c("No English")

constraint_list[[6]] <- list()
constraint_list[[6]][["Faces"]] <- c("White man", "White woman", 
                                     "Light man", "Light woman",
                                     "Medium man", "Medium woman")
constraint_list[[6]][["Identifies"]] <- c("Black")

constraint_list[[7]] <- list()
constraint_list[[7]][["Faces"]] <- c("Brown man", "Brown woman",
                                     "Black man", "Black woman")
constraint_list[[7]][["Identifies"]] <- c("White")

marginal_weights <- list()
marginal_weights[["Faces"]] <- rep(1/length(attribute_list[["Faces"]]),
                                   length(attribute_list[["Faces"]]))
marginal_weights[["Background"]] <- rep(1/length(attribute_list[["Background"]]),
                                        length(attribute_list[["Background"]]))
marginal_weights[["Identifies"]] <- rep(1/length(attribute_list[["Identifies"]]),
                                        length(attribute_list[["Identifies"]]))
marginal_weights[["Language"]] <- rep(1/length(attribute_list[["Language"]]),
                                      length(attribute_list[["Language"]]))
marginal_weights[["Education"]] <- rep(1/length(attribute_list[["Education"]]),
                                       length(attribute_list[["Education"]]))
marginal_weights[["Profession"]] <- rep(1/length(attribute_list[["Profession"]]),
                                        length(attribute_list[["Profession"]]))
marginal_weights[["Religion"]] <- rep(1/length(attribute_list[["Religion"]]),
                                      length(attribute_list[["Religion"]]))
marginal_weights[["Political"]] <- c(2/9, 2/9, 2/9, 3/9)
marginal_weights[["Spouse"]] <- c(5/15, 3/15, 2/15, 2/15, 3/15)

conjointdesign <- cjoint::makeDesign(type='constraints', attribute.levels=attribute_list,
                                     constraints=constraint_list, level.probs=marginal_weights)
baselines <- list()
baselines$Education <- "High school graduate"
baselines$Background <- "US-born legal grandparents"
baselines$Profession <- "Teacher"
baselines$Political <- "Not political"
baselines$Religion <- "Not religious"

####Feature order -----
c_order <- c("Faces", "Identifies", "Religion", "Background", "Language",
             "Spouse", "Education", "Profession", "Political") 

#### AMCE of Perceptions of Americanness (Amerian rating) ----
####NOTE: This model is computationally intensive and may not run on all computers
ratingamce_r <- cjoint::amce(rating_r ~ Faces + Background + Identifies +
                               Language + Education + Profession + Religion 
                             + Political + Spouse, data=df, cluster=TRUE, 
                             respondent.id="caseid", design=conjointdesign,
                             weights="weight", baselines=baselines,
                             na.ignore=TRUE)


#### Robustness for AMCE perceptions of americanness among all respondents ----
amce.sum <- cjoint::summary.amce(ratingamce_r) ## extract values

## Uniform mixture prior 
beta_ash <- amce.sum$amce[3][,1] ## get beta.hat
se.beta_ash <- amce.sum$amce[4][,1] ## get se.hat

ashR.unif <- ash(beta_ash, se.beta_ash, ## supply the beta.hat & se.hat
                 mixcompdist= "uniform", ## specify the spike-slap prior, 
                 ## "uniform" mixture here
                 method = "fdr")$result  ## "fdr" b/c main goal is to assess
## false discovery rate

ashR.unif.post.beta <- ashR.unif[,9] ## get the posterior beta.hat
ashR.unif.post.sd <- ashR.unif[,10] ## get the posterior se.hat
ashR.unif.post.qVal <- ashR.unif[,8] ## get the q value from false disc. rate

## Normal mixture prior 

ashR.norm <- ash(beta_ash, se.beta_ash, 
                 mixcompdist = "normal", ## similar to above, but "normal" 
                 ## mixture here
                 method = "fdr")$result

ashR.norm.post.beta <- ashR.norm[,9]
ashR.norm.post.sd <- ashR.norm[,10]
ashR.norm.post.qVal <- ashR.norm[,8]


ashR.norm.post.beta

###Table A16----
out_amce_robust <- tibble(attributes = amce.sum$amce$Attribute,
       levels = amce.sum$amce$Level,
       new.b = ashR.norm.post.beta, 
       new.sd = ashR.norm.post.sd, 
       new.q = ashR.norm.post.qVal,
       new.star = case_when(new.q < 0.01 ~ "**", 
                            new.q < 0.05 ~ "*", 
                            new.q < 0.1 ~ "+", 
                            T ~ ""),
       amce.b = amce.sum$amce$Estimate, 
       amce.sd = amce.sum$amce$`Std. Err`, 
       amce.p = amce.sum$amce$`Pr(>|z|)`, 
       amce.star = case_when(amce.p < 0.01 ~ "**", 
                             amce.p < 0.05 ~ "*", 
                             amce.p < 0.1 ~ "+", 
                             T ~ "")) %>% 
  select(-contains(c(".p", ".q")))

print.xtable(xtable(out_amce_robust), "d_tables/table_A16.tex", type = "latex", include.rownames=FALSE)


#### MM diff of perceived americanness, between W and B respondents ----

american_diff <- cregg::mm_diffs(na.omit(df),
                                 rating_r ~ Faces + Background
                                 + Education + Profession
                                 + Language + Political
                                 + Religion + Spouse
                                 + Identifies,
                                 by = ~ race,
                                 weights = ~weight, 
                                 feature_order = c_order)

## Uniform mixture prior 
beta_ash_b <- american_diff$estimate ## get beta.hat
se.beta_ash_b <- american_diff$std.error ## get se.hat

ashR.unif_b <- ash(beta_ash_b, se.beta_ash_b, ## supply the beta.hat & se.hat
                 mixcompdist= "uniform", ## specify the spike-slap prior, 
                 ## "uniform" mixture here
                 method = "fdr")$result  ## "fdr" b/c main goal is to assess
## false discovery rate

ashR.unif.post.beta_b <- ashR.unif_b[,9] ## get the posterior beta.hat
ashR.unif.post.sd_b <- ashR.unif_b[,10] ## get the posteria se.hat
ashR.unif.post.qVal_b <- ashR.unif_b[,8] ## get the q value from false disc. rate

## Normal mixture prior 

ashR.norm_b <- ash(beta_ash_b, se.beta_ash_b, 
                 mixcompdist = "normal", ## similar to above, but "normal" 
                 ## mixture here
                 method = "fdr")$result

ashR.norm.post.beta_b <- ashR.norm_b[,9]
ashR.norm.post.sd_b <- ashR.norm_b[,10]
ashR.norm.post.qVal_b <- ashR.norm_b[,8]


ashR.norm.post.beta_b

###Table A17----
out_americanbyrace_robust<- tibble(attributes = american_diff$level,
                          levels = american_diff$Level,
                          new.b = ashR.norm.post.beta_b, 
                          new.sd = ashR.norm.post.sd_b, 
                          new.q = ashR.norm.post.qVal_b,
                          new.star = case_when(new.q < 0.01 ~ "**", 
                                               new.q < 0.05 ~ "*", 
                                               new.q < 0.1 ~ "+", 
                                               T ~ ""),
                          mmdiff.b = american_diff$estimate, 
                          mmdiff.sd = american_diff$std.error, 
                          mmdiff.p = american_diff$p, 
                          mmdiff.star = case_when(mmdiff.p < 0.01 ~ "**", 
                                                  mmdiff.p < 0.05 ~ "*", 
                                                  mmdiff.p < 0.1 ~ "+", 
                                                T ~ "")) %>% 
  select(-contains(c(".p", ".q")))

print.xtable(xtable(out_americanbyrace_robust), "d_tables/table_A17.tex", type = "latex", include.rownames=FALSE)


#### MM diff by American identity strength among White respondents ----

amerimport_diff_w <- mm_diffs(na.omit(white_long),
                              rating_r ~ Faces + Background
                              + Education + Profession
                              + Language + Political
                              + Religion + Spouse
                              + Identifies,
                              id = ~ caseid,
                              by = ~ amerimport_high,
                              weights = ~weight,
                              feature_order = c_order)

#### Robustness check for marginal mean differences, by American identity strength, among W respondents ----

## Uniform mixture prior 
beta_ash_c <- amerimport_diff_w$estimate ## get beta.hat
se.beta_ash_c <- amerimport_diff_w$std.error ## get se.hat

ashR.unif_c <- ash(beta_ash_c, se.beta_ash_c, ## supply the beta.hat & se.hat
                   mixcompdist= "uniform", ## specify the spike-slap prior, 
                   ## "uniform" mixture here
                   method = "fdr")$result  ## "fdr" b/c main goal is to assess
## false discovery rate

ashR.unif.post.beta_c <- ashR.unif_c[,9] ## get the posterior beta.hat
ashR.unif.post.sd_c <- ashR.unif_c[,10] ## get the posteria se.hat
ashR.unif.post.qVal_c <- ashR.unif_c[,8] ## get the q value from false disc. rate

## Normal mixture prior 

ashR.norm_c <- ash(beta_ash_c, se.beta_ash_c, 
                   mixcompdist = "normal", ## similar to above, but "normal" 
                   ## mixture here
                   method = "fdr")$result

ashR.norm.post.beta_c <- ashR.norm_c[,9]
ashR.norm.post.sd_c <- ashR.norm_c[,10]
ashR.norm.post.qVal_c <- ashR.norm_c[,8]


ashR.norm.post.beta_c

####Table A18----
out_americanidwhite_robust<- tibble(attributes = amerimport_diff_w$level,
                                   levels = amerimport_diff_w$Level,
                                   new.b = ashR.norm.post.beta_c, 
                                   new.sd = ashR.norm.post.sd_c, 
                                   new.q = ashR.norm.post.qVal_c,
                                   new.star = case_when(new.q < 0.01 ~ "**", 
                                                        new.q < 0.05 ~ "*", 
                                                        new.q < 0.1 ~ "+", 
                                                        T ~ ""),
                                   mmdiff.b = amerimport_diff_w$estimate, 
                                   mmdiff.sd = amerimport_diff_w$std.error, 
                                   mmdiff.p = amerimport_diff_w$p, 
                                   mmdiff.star = case_when(mmdiff.p < 0.01 ~ "**", 
                                                           mmdiff.p < 0.05 ~ "*", 
                                                           mmdiff.p < 0.1 ~ "+", 
                                                           T ~ "")) %>% 
  select(-contains(c(".p", ".q")))
print.xtable(xtable(out_americanidwhite_robust), "d_tables/table_A18.tex", type = "latex", include.rownames=FALSE)


#### MM diff by racial identity strength among Black respondents ----

raceimport_diff_b <- mm_diffs(na.omit(black_long),
                              rating_r ~ Faces + Background
                              + Education + Profession
                              + Language + Political
                              + Religion + Spouse
                              + Identifies,
                              id = ~ caseid,
                              by = ~ raceimport_high,
                              weights = ~weight, 
                              feature_order= c_order)

#### Robustness check for marginal mean differences, by racial identity strength, among B respondents ----

## Uniform mixture prior 
beta_ash_d <- raceimport_diff_b$estimate ## get beta.hat
se.beta_ash_d <- raceimport_diff_b$std.error ## get se.hat

ashR.unif_d <- ash(beta_ash_d, se.beta_ash_d, ## supply the beta.hat & se.hat
                   mixcompdist= "uniform", ## specify the spike-slap prior, 
                   ## "uniform" mixture here
                   method = "fdr")$result  ## "fdr" b/c main goal is to assess
## false discovery rate

ashR.unif.post.beta_d <- ashR.unif_d[,9] ## get the posterior beta.hat
ashR.unif.post.sd_d <- ashR.unif_d[,10] ## get the posteria se.hat
ashR.unif.post.qVal_d <- ashR.unif_d[,8] ## get the q value from false disc. rate

## Normal mixture prior 

ashR.norm_d <- ash(beta_ash_d, se.beta_ash_d, 
                   mixcompdist = "normal", ## similar to above, but "normal" 
                   ## mixture here
                   method = "fdr")$result

ashR.norm.post.beta_d <- ashR.norm_d[,9]
ashR.norm.post.sd_d <- ashR.norm_d[,10]
ashR.norm.post.qVal_d <- ashR.norm_d[,8]


ashR.norm.post.beta_d

###Table A19----
out_raceidblackrobust<- tibble(attributes = raceimport_diff_b$level,
                               levels = raceimport_diff_b$Level,
                               new.b = ashR.norm.post.beta_d, 
                               new.sd = ashR.norm.post.sd_d, 
                               new.q = ashR.norm.post.qVal_d,
                               new.star = case_when(new.q < 0.01 ~ "**", 
                                                    new.q < 0.05 ~ "*", 
                                                    new.q < 0.1 ~ "+", 
                                                    T ~ ""),
                               mmdiff.b = raceimport_diff_b$estimate, 
                               mmdiff.sd = raceimport_diff_b$std.error, 
                               mmdiff.p = raceimport_diff_b$p, 
                               mmdiff.star = case_when(mmdiff.p < 0.01 ~ "**", 
                                                       mmdiff.p < 0.05 ~ "*", 
                                                       mmdiff.p < 0.1 ~ "+", 
                                                       T ~ "")) %>% 
  select(-contains(c(".p", ".q")))

print.xtable(xtable(out_raceidblackrobust), "d_tables/table_A19.tex", type = "latex", include.rownames=FALSE)


#### Clear environment ----

rm(list = ls())

